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r~| \ Abstract. Anomalies in direct and indirect detection have motivated models of dark matter 

consisting of a multiplet of nearly-degenerate states, coupled by a new GeV-scale interaction. 
, We perform a careful analysis of the thermal freezeout of dark matter annihilation in such 

a scenario. We compute the range of "boost factors" arising from Sommerfeld enhancement 
in the local halo for models which produce the correct relic density, and show the effect of 
including constraints on the saturated enhancement from the cosmic microwave background 
(CMB). We find that boost factors from Sommerfeld enhancement of up to ~ 800 are possible 
in the local halo. When the CMB bounds on the saturated enhancement are applied, the 
maximal boost factor is reduced to ~ 400 for 1-2 TeV dark matter and sub-GeV force 
carriers, but remains large enough to explain the observed Fermi and PAMELA electronic 
signals. We describe regions in the DM mass-boost factor plane where the cosmic ray data is 
well fit for a range of final states, and show that Sommerfeld enhancement alone is enough 
to provide the large annihilation cross sections required to fit the data, although for light 
mediator masses (ms < 200 MeV) there is tension with the CMB constraints in the absence 
of astrophysical boost factors from substructure. Additionally, we consider the circumstances 
under which WIMPonium formation is relevant and find for heavy WIMPs (~ 2 TeV) and 
soft-spectrum annihilation channels it can be an important consideration; we find regions 
with m x > 2.8 TeV that are consistent with the CMB bounds with 0(600-700) present-day 
boost factors. 
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1 Introduction 

Recent measurements of electron and positron cosmic rays at 10-1000 GeV [1-7] indicate a 
rise in the positron flux fraction at 10-100+ GeV and a hardening of the e + + e~ spectrum at 
~ 20— 1000 GeV, suggesting a new source of positrons and electrons with a hard spectrum and 
a TeV-scale cutoff. Annihilation of weak-scale dark matter provides a natural explanation for 
an excess of pairs at this energy scale, but conventional models of thermal relic dark matter 
annihilation are challenged by the large amplitude and hard spectrum of the signal [8, 9], as 
well as the absence of any corresponding excess in cosmic ray antiprotons [9, 10]. 

The presence of a new GeV-scale force in the dark sector has been proposed to explain 
these features of the measured cosmic-ray excesses [11-13]. The dark matter can efficiently 
annihilate to the force carriers (f), which then decay into pairs of Standard Model particles; 
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if < 2m p , the decay into antiprotons is kinematically forbidden, and the resulting e + 
spectrum is hard due to the highly boosted intermediate (on-shell) light state [14, 15] . Addi- 
tionally, 4> mediates a fm-range attractive DM self-interaction, which enhances annihilation 
at low velocities, potentially by several orders of magnitude, i.e., the Sommerfeld enhance- 
ment [16] . The presence of 0(MeV) splittings between dark matter states can be used to 
explain the 511 keV excess in the Galactic center observed by INTEGRAL/SPI [19, 20], via 
collisional excitation and decay of the excited state ("eXciting dark matter" (XDM)) [14, 21- 
24], and to reconcile the annual modulation observed by DAMA/LIBRA [25-27] with the 
null results of other direct detection experiments through inelastic WIMP-nuclear scattering 
("inelastic dark matter" (iDM)) [28, 29]. 

While the low-velocity enhancement to annihilation has been well explored, and fits to 
the existing data from the new annihilation channels have been given [30, 31], no complete 
picture has yet been presented. The Sommerfeld enhancement is strongly velocity-dependent 
(scaling as 1/v generically) , and thus has a much larger effect in the local halo than during 
freezeout, when the dark matter thermal relic density is established. However, the effect on 
freezeout can still be significant: an 0(1) enhancement to the annihilation rate at freezeout 
leads to 0(1) changes in the thermal relic density, so the underlying "bare" annihilation 
cross section must be reduced to compensate [32, 33]. In narrow regions of parameter space, 
where the enhancement scales as 1/v 2 rather than 1/v below a certain velocity, this effect 
can be more pronounced: in particular, in such "resonance regions" it may be possible for 
DM annihilation to "un-freeze" after kinetic decoupling (once the DM begins to cool faster 
than the CMB), leading to a very efficient depletion of the relic density [11, 34]. Indeed, in 
cases where XX ~~ 4> ~^ there is significant tension between these bounds and the 

annihilation rates required to generate the cosmic-ray anomalies [34, 35]. 

In this work, we study a broader range of parameter space, including a range of decay 
modes for the <p and splittings between dark matter states. We will see that while the 
scenarios discussed by [34] are tightly constrained, generic models have less tension due 
to enhanced annihilation in the presence of mass splittings and the typical hard electron 
component arising from <p — > e + e~. We will employ the approximations derived in [36] 
to identify regions of interest and present spectra for self-consistent benchmark points that 
fit the cosmic-ray data; we verify that a careful numerical calculation of the Sommerfeld 
enhancement at these benchmark points does not appreciably alter our results. We find 
that the most stringent constraints on the present-day annihilation cross section originate 
not from achieving the appropriate relic abundance, but from WMAP measurements of the 
CMB temperature and polarization angular power spectra [37, 38], which will be greatly 
improved by Planck. We motivate and describe the models we consider in §2, outline the 
details of our calculation in §3, and discuss constraints from indirect detection in §4. In §5 
we describe our modeling of the cosmic ray spectra, and our criterion for a "good fit" to the 
cosmic ray data. Readers principally interested in our results can immediately skip to §6, 
which is essentially self-contained. 

Throughout this work we will use the term "boost factor" (BF) to mean the s-wave 
annihilation cross section {av} divided by 3 x 10~ 26 cm 3 /s (this value has been employed in 
the literature as the canonical value of {av} for thermal relic DM). 



1 This effect was originally studied in the context of dark matter annihilation in [17, 18]. 
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2 Inelastic Dark Matter Through the Vector Portal 

Let us begin by laying out the model space we consider. We focus here on models where 
dark matter freezes out by annihilating into a new force carrier XX ~~ 'P'P [14, 39], which 
itself maintains kinetic equilibrium with the standard model, typically through the process 
e7 -H- e<j) [40], or s-channel dark boson exchange [41]. While models can be constructed 
through a variety of "portals", we focus our attention on the vector portal. That is, we 
assume <fi is the vector boson of a new gauge group U(1)d, Higgsed at roughly the GeV scale, 
and the connection between dark forces and the standard model comes through an effective 
term eFj^^Fj^ . The dark sector is generically thermalized before dark matter freezeout for 

e ~ 10~ 7 , although it should be noted that so long as the sectors are brought into thermal 
contact through some other physics earlier, this condition could be relaxed [39]. 

It is important to emphasize that because we are considering a vector force, we are 
forced to consider sectors with multiple states. The minimal fermionic representation of 
U(1)d is a Dirac fermion, which is composed of two Majorana fermions, while the minimal 
scalar representation is complex, composed of two real scalars. Once U(1)d is broken, absent 
an accidental low-energy global symmetry, there is no a priori reason these states should 
remain degenerate. If the components are degenerate, then direct detection experiments 
[42, 43] constrain e ~ 10~ 6 (10~ 8 ) for rru = 1 GeV (100 MeV). In supersymmetric models, 
oc e, yielding a fairly robust cross section of ~ 10~ 38 cm 2 [44], well above current direct 
detection limits. However, if the different WIMP components are non-degenerate by an 
amount 5 ~ m x v 2 , then because the vector coupling is off-diagonal [45], the elastic scattering 
cross section will be suppressed or eliminated. 

Thus we should emphasize that within this overall framework the most generic setup 
is one where dark matter consists of multiple states (i.e., a pseudo-Dirac fermion), split by 

an amount ~ 100 keV, and where the force carrier interacts with SM matter through its 
mixing with the photon. This final point is especially important, because as we shall see, 
the splitting naturally enhances the late-time Sommerfeld enhancement, and the coupling 
of 4> to charge generally produces a sizeable (ft — > e + e~ component, which is the dominant 
contributor to the DM explanation of signals seen at PAMELA and Fermi. Other models 
can be constructed, for instance through the Higgs portal or through the axion portal [13], 
but degenerate WIMPs annihilating dominantly into 4// is actually extremely challenging in 
the context of vector portal models, and not generically realized within this setup. 

2.1 Annihilation channels and Sommerfeld enhancement 

As our canonical example, we consider a model where the dark matter is a Dirac fermion 
charged under a hidden sector U(l)n, broken at the GeV scale by an Abelian Higgs. At low 
energies, we assume higher dimension operators will split the Majorana components xi,2 of 
the WIMP by an amount m X2 — m Xl = 5 ~ MeV. Even at high energies, it is convenient to 
think of annihilation processes in terms of these states. 

Most generically, the dominant annihilation is XiXi> X2X2 —> <fi<l>, via t-channel exchange 
of the other x state. There are also coannihilation channels into Higgses (which, at late times 
can be cast in the unitary gauge as X1X2 — > 4>h-D), or other light states charged under U(1)d- 
Both of these channels experience Sommerfeld enhancement. The unenhanced tree- level 
annihilation cross section into two gauge bosons, for nonrelativistic or mildly relativistic 
fermions, is given by the usual result for pair production (e.g. [46]), with the dominant low- 
velocity contribution coming from the s-wave term, cx| f rc i | z=o = na 2 D /m x . Here v re \ is the 
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relative velocity of the two interacting particles, and p is the 3- momentum of a single initial 
fermion in the COM frame. 

The cross section for s-channel coannihilation into two charge- 1 dark Higgses is given 

by, 



\p\ (m 2 , + \p\ 



relU=0 



and scales as the Higgs charge squared. The leading s-wave contribution is a\v } 
trajj 1 4m^. , a factor of 4 lower than for the ^-channel annihilation (for a singly charged Higgs) . 

To calculate the Sommerfeld enhancement for these channels, one must sum a series of 
ladder diagrams (or equivalently, solve the Schrodinger equation). The interaction between 
the dark gauge bosons and the Xi eigenstates is off-diagonal, so an initial two-particle state 
XiXi can be scattered into the two-particle state X2X2 (or vice versa), but not into the state 
XiX2- Similarly, for particles initially in a two-body X1X2 state, the only effect of the long- 
range vector interaction is to swap the individual particle states. Thus, when computing the 
potential due to the long-range interaction, the X1X2 2-body state is disjoint from the other 
two and can be treated separately. 

Since the vector-mediated scatterings are always purely elastic for the X1X2 2-body 
state, and the long-range interaction is attractive, the effect of the Sommerfeld enhancement 
is well approximated by the enhancement due to a Yukawa potential (up to corrections of 
order 5 /m x , due to the slightly different masses of the interacting particles) . 

For the X1X1 an d X2X2 states, the interaction at tree-level is always inelastic. The 
corresponding scattering problem can be written in terms of a potential matrix with off- 
diagonal Yukawa terms, and approximately solved as in [36]. We employ the approximate 
semi-analytic results of [36] to estimate the Sommerfeld enhancement; for our benchmark 
points, we also check the approximation numerically. 

2.2 A specific model for the mass splitting 

In the case where the Higgs is singly charged, an O(MeV) mass splitting can be generated 
naturally by a higher-dimension operator, which in the high-energy theory (where the U(1)d 
is unbroken) takes the form (l/2)(y/A)(^ c *$>h* D h* D + h.c). In Appendix B we calculate the 
mass splitting and annihilation rate arising from this operator; here we will summarize those 
results. 

In terms of the xi,2 mass eigenstates, the operator takes the form, 

Aiplit = (!/4)(y/A)(xiXl(^I)^D + /iD/iD)-X2X2(/iD/lD + /iD/iD)+2iXlX2(/lD/iD-/iD/iD))+/l-C. 

Working in unitarity gauge for simplicity, and writing hp — > (vd + /id)/v / 2, we obtain, 
Aplit -> (l/4)(y/A)(xiXi - X2X2){v 2 D + h 2 D + 2v D h D ) + h.c. 

The Yukawa coupling is suppressed by the small ratio vd/A, but the (l/4)(y/A)xiXi^i> term 
induces a potentially large XiXi annihilation channel, as discussed in [47]. Furthermore, the 
mass splitting 5 = {y/K)v 2 D , and the mediator mass = g£> V£>: this allows us to write the 
low-velocity cross section for XiXi ~^ ^d^d annihilation as, 

1 9 / <5m v \ ira 2 n 

= 2" {-4 ) -4- < 2 - 2 > 
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This relation comes from a tree-level computation and neglects Sommerfeld corrections. 
Since the XiXi ^d^d channel corresponds to annihilation of two same-charge fermions in 
the high-energy limit, it is actually Sommerfeld suppressed: in the ladder diagram picture, 
this can be understood as due to destructive interference between the two relevant diagrams 
(with the annihilation being through X1X1 or X2X2 in the final step of the ladder diagram), 
originating from the fact that the XiXi^d^d and X2X2 Tid^d couplings have opposite signs. 
Since Sommerfeld-suppressed channels are negligible for low v, we will approximate the Som- 
merfeld suppression by the result for p-wave Coulomb scattering, 

Srep.p-wave = (^/l _ f 1 + 4^2) ' ( 2 ' 3 ) 

This channel is therefore always negligible in the present day, as it experiences both a 
Sommerfeld suppression and a p- wave suppression. However, it can be important at freezeout, 
depending on the ratio m x 5/m^. If the DM were instead a complex scalar, the p-wave 
suppression would be absent and the effect on freezeout much larger 2 . 

3 Solving for the Relic Density 

We solve numerically for the abundances in the ground (1) and excited (2) states, using the 
publicly available IDL code LSODE. We define Y{ = rii/s, where s is the entropy density, and 
x = rn-y/T. The Boltzmann equation becomes, 

it = Htk) {FY2 " S " Y ^ {a ^ v) + ( ™ ~ Y ^^ + ~ , 

dY 2 _ X , r/ v 2 v2w„A„,\ , / V .V v2w„A„,\ , u_ v 2 



dx H(m y 



-TY 2 - s [(Y 2 2 - Y c l){a*v) + (Y,Y 2 - Y c \){a^v) + k D Y 2 2 - k E Y?}) . (3.1) 



Here T is the decay rate for the excited state X2> and k E and ku describe (respectively) 
the excitation and de-excitation of the excited state by DM-DM scattering. We include 
the s- and p-wave contributions to the annihilation cross sections. We have verified that 
including the complete (tree-level) relativistic cross sections, with the s- and p-wave pieces 
corrected by Sommerfeld enhancement, and using the full relativistic momentum distribution 
for thermal particles (rather than the nonrelativistic Maxwell-Boltzmann distribution), does 
not significantly alter our results. 

y^ q , the equilibrium value of the abundances (normalized to the entropy density), is 
determined simply by the number of degrees of freedom of the DM compared to the SM 
relativistic degrees of freedom, and the temperature and mass of the DM. Strictly the ground 
and excited states have distinct values of Y eq , but due to the large hierarchy between the 
dark symmetry breaking scale and the dark matter mass, the splitting between the states 
is irrelevant during freezeout. By the time the temperature of the universe drops to the 
symmetry breaking scale, much less the scale of the mass splittings, is infinitesimal and 
irrelevant to the evolution of Y\, Y 2 . Thus it is acceptable to approximate Y^ q = Y^. 

The Weyl fermions xi an d xi each have g = 2 internal degrees of freedom, and we use 
the standard results 



JLf^JL-x^e-*, x»3, 

x < 3. 



We thank Josh Ruderman for this observation. 



- 5 - 



Here g e s = 3g/4 for the fermionic dark matter we consider; for bosonic DM g c g = g. g*s 
counts the relativistic degrees of freedom for the entropy density. 

Once the DM thermally decouples from the SM photon bath, it cools more rapidly 
and the Sommerfeld enhancement becomes more pronounced [11, 32, 33]. Prior to kinetic 

— 1 /3 

decoupling, the DM temperature is equal to the CMB temperature, scaling as g^ s /a; after 
decoupling, it scales as 1/a 2 . Thus we obtain the relation, 



Tdm 



T CMB ( 9*s{TcMB] 



Trd V 9*s{Trd) 



2/3 



We employ the expression for the kinetic decoupling temperature derived in [34], with 
the added requirement that the kinetic decoupling temperature must be larger than the mass 
splitting, so that inelastic scatterings of DM on SM particles are not kinematically suppressed, 



- kd 



~ max < m. 



5, 0.82MeV 



10' 



1/2 



30MeV 



0.021 

OLD 



1/4 



m X 

TeV 



1/4 



(3.3) 



We take e = 10 -3 as a benchmark, but our results are not very sensitive to this choice. 

We have neglected the terms in the Boltzmann equation describing the up- and down- 
scattering of DM states on SM particles: at temperatures where the mass splitting becomes 
significant and the ground and excited states may have different populations, these processes 
are (1) slow compared to a Hubble time, and (2) subdominant compared to the effects of 
DM-DM interactions, in models with Sommerfeld-enhanced scattering. (1) follows simply 
from the observation that kinetic decoupling occurs at temperatures higher than the mass 
splitting; (2) can be seen e.g. from [34], where DM-DM interactions are found to become 
inefficient at preserving thermal equilibrium at ~ 10 keV, well below the kinetic decoupling 
temperature. 

The excitation fraction affects the annihilation rate in two ways. The smaller effect is 
that the Sommerfeld enhancements differ for XiXi> X2X2 and X1X2 initial states (although 
of course they are identical in the limit where the temperature greatly exceeds the mass 
splitting). 

More importantly, the annihilation cross section for a X1X2 initial state differs from 
that for a X1X1 or X2X2 initial state, so the annihilation cross section averaged over all 
possible pairs of interacting DM particles is a function of the excitation fraction. Let the 
total DM number density be denoted n, and the number densities for the ground and excited 
states be denoted n\ and n<i respectively. Then the number of annihilations per volume 
per second, in the absence of Sommerfeld enhancement, is given for this simple model by 
n\(pY]V) /2 + n^^a^y) /2 + n\nii<J\i v ) ■ The first two terms come from t-channel annihilation 
into 0's, the last term from s-channel annihilation into dark Higgses. Let us write a = 
Ct\ = = ^&\2-i taking the s-wave annihilation cross sections derived previously. Then at 



'11 — "22 
freezeout, n\ 



77-2 ~ n/2 and the total annihilation rate per unit volume is given by, 



(H((n/2) 2 + (l/4)(n/2) 2 ) = -Mn 2 . 

In the present day, 122 ~ and n\ ~ n, so the annihilation rate per unit volume is (av)n 2 /2] 
larger by a factor of k = 8/5. More generally, a lower excitation fraction leads to a greater 
annihilation rate. However, this statement is entirely a function of the ratio of the xiXij 
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X1X2 and X2X2 annihilation cross sections: in particular, if there are other states present 
in the dark sector that couple to the dark gauge boson, the s-channel X1X2 annihilation 
is significantly enhanced due to the greater number of final states, leading to a reduction 
or even reversal of this effect. If the s-channel X1X2 annihilation cross section is given by 
( a i2 v )i=o = Kirajj/m^., the rescaling factor k becomes k = 2/(1 + k). Note, however, that 
such additional or enhanced X1X2 annihilation channels also modify the value of required 
to obtain the correct relic density, and therefore change the Sommerfeld enhancement: the 
present-day annihilation cross section cannot simply be rescaled by k relative to the case 
with k = 1. We will focus on the minimal k = 1/4 case corresponding to a singly charged 
Higgs, with a xxhh operator generating the mass splitting (and giving rise to an additional 
annihilation channel) as described in §2.2, but we will also present results for k = 1 and 4 
for illustrative purposes. For these latter cases we do not introduce an explicit mechanism to 
generate the mass splittings; we assume the only contribution to the self-annihilation cross 
section arises from the usual XX ~ > 4'4> channel, and parameterize the coannihilation cross 
section by k. 

3.1 DM-DM scattering 

DM-DM scattering does not affect the relic density of dark matter directly, but when the 
temperature of the universe drops below the mass splitting, the excitation and de-excitation 
rates can affect the relative populations of the ground and excited states, which in turn 
impacts the annihilation rate. However, because of the large hierarchy between the DM mass 
and the mass splittings, these effects only come into play long after freezeout, and thus have 
little effect on the relic density except very close to resonances. However, we include them 
for completeness. 

We are interested in regions of parameter space where a large Sommerfeld enhancement 
persists at low velocities; consequently, the mass splitting cannot be much greater than 
the Rydberg energy of the XX system [36]. It follows that at low velocities the dominant 
contribution to the DM-DM scattering arises from infinite ladder diagrams; the principal 
differences between elastic and inelastic scattering (or excitation vs de-excitation) are then 
just the transferred momentum (q 2 > m x S for inelastic scattering) and the final-state phase 
space (since the ladders differ only in the masses of the initial- and final-state particles, and 
we assume 5 <C m x ). Therefore we expect the elastic scattering, excitation and de-excitation 
cross sections to be roughly equal, up to phase space factors and with the replacement 
r?i0 — > y/m x 5 if m x v < y/m x 5 < m^. For the elastic scattering cross section we employ the 
prescription of [34] for the momentum transfer cross section, based on [49], 



where /3 = 2aD'm ( j > /(m x v 2 el ). 

3.2 Decay of the excited state 

In our current toy model, for mass splittings much smaller than 2m e , the lifetime of the 
excited state exceeds the age of the universe [24]. This is potentially problematic due to 
constraints from direct detection on the present-day excitation fraction (since the excited 




P<0.1, 
'<), 0.1 < p < 1000, 

/3) 2 , /3> 1000, 



(3.4) 



7T 



(ln/3 + 1 
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state can downscatter in detectors) [24]: consequently, for 5 <C 2m e , either the excited 
state must be efficiently depleted by DM-DM scattering, or the model must have some extra 
ingredient that can mediate the decay of the excited state. A lifetime between 10 13 — 10 18 
s may be of particular interest, as that would allow the annihilation rate at the redshift 
of last scattering to be reduced relative to the present day, since the DM annihilates more 
efficiently when it is entirely in the ground state (at least for models with k < 1), weakening 
the CMB constraints. Once the mass splitting exceeds 2m e , opening up decay to an e + e~ 
pair, the lifetime of the excited state is naturally around that of the neutron, or ~ 10 3 s for 
e ~ I0~ 3 (the lifetime scales as 1/e 2 ) [24]. In any case, these lifetimes greatly exceed the age 
of the universe at kinetic decoupling, and cannot significantly affect the relic density. For 
the purposes of this calculation, we set T = 6.6 x 10 -28 GeV, corresponding to a lifetime of 
10 3 s, as appropriate for a mass splitting greater than 2m e with e ~ 10 -3 . We have verified 
that our results are insensitive to the decay lifetime, at least for the relatively short lifetimes 
relevant for 5 > 2m e : increasing the lifetime up to 10 9 s (corresponding to e ~ 10~ 6 ) has no 
measurable impact on our results, and increasing the lifetime to infinity affects our results 
only in the sense that annihilation in the present day can also involve excited-state DM (that 
is, it affects the present-day annihilation rates relevant for indirect detection, but does not 
modify the relic density). 

4 Constraints from the Cosmic Microwave Background and Other Indirect 
Searches 

Many bounds on dark matter annihilation from indirect detection searches have been dis- 
cussed in the literature. However, most of them have large uncertainties associated with 
model-dependent or poorly known astrophysical factors, and so we do not employ them as 
constraints on our allowed regions of parameter space; later in this section we will briefly out- 
line these bounds and justify neglecting them. The exception is a set of constraints from the 
cosmic microwave background that provide uniquely clean bounds on the dark matter anni- 
hilation cross section in models with Sommerfeld-enhanced annihilation, with no dependence 
on Galactic astrophysics or the history of DM structure formation, and minimal dependence 
on the details of the DM model; we now give a brief outline of these constraints and their 
application. 

Dark matter annihilation around the redshift of last scattering can significantly perturb 
the temperature and polarization angular power spectra of the cosmic microwave background 
(CMB) [37, 38, 50] 3 . DM models with Sommerfeld-enhanced annihilation are especially 
sensitive to the resulting constraints, since the annihilation cross section increases at low 
velocity. The 95% confidence limits on DM annihilation from WMAP5 can be expressed as 
[38] 4 , 

jj^oM < 120 / m x \ 
3 x 10- 26 cm 3 /s ~ / VlTeV/ ' 1 ' ' 

Here / is a parameter describing the fraction of energy from DM annihilation which ionizes 
and heats the intergalactic medium; [38] showed that / can be approximated as / ~ 0.7 for 
annihilation to electrons and / ~ 0.2 — 0.3 for all other SM final states, excepting neutrinos, 

3 Limits from the CMB spectrum itself are much weaker than the temperature and polarization anisotropy 
constraints [33, 51]. 

4 Note that these constraints assume a constant primordial scalar spectral index n s in re-fitting the cosmo- 
logical parameters; if running of the spectral index is allowed, these constraints may be weakened somewhat. 
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where / ~ (that is, / ~ if the DM annihilates directly to neutrinos; if the DM annihilates 
to unstable SM particles which then decay to neutrinos, the value / ~ 0.2 — 0.3 should 
be used). Thus provided DM annihilations have no significant branching ratio directly to 
neutrinos, the boost factor at very low velocity is constrained to satisfy BF(u — > 0) < 
600?n x /l TeV, with a smaller limit if the dark sector annihilation products decay to electrons 
with a non-negligible branching ratio. 

This constraint can be combined with the semi-analytic expression for the Sommerfeld 
enhancement to obtain a DM-mass-independent limit on the allowed boost factor in the 
present day. Using the semi-analytic form for the Sommerfeld enhancement derived in [36], 
as v — > the enhancement saturates at, 

~ \ fi J 1 - cos {e 5 7r/n + 26L) > fi ' ' ' ' 



Here /x is a function of = m^j ctDm x and e$ = x/25/m x /aD (defined in [36]), but in 
practice it is nearly independent of e^: for the region where the approximations of [36] are 
expected to be accurate (5 <C ar/m^ <C a 2 D m x ), we find fi « (1/2)(1 + \fh)m^j au^x t° 
within 30%. 

In the present day, if fi <C e^, then the enhancement is bounded above by S < 2n/e v (in 
the inelastic case, S ~ 2ir/e v , whereas in the elastic ir/e v ). Then the ratio of the 

present-day enhancement to the saturated enhancement is bounded by, 



Snow ^ 2lT I € v 2 [ M \ ^ 



'sat 



7T 2 //i TT \€ v J TT 



(4.3) 



where = ar>m x \i. The limit from WMAP5 requires that BF sat < (120//)m x /lTeV, 
which in turn yields a limit on the present-day boost factor depending only on m<*, 5 and the 
local velocity dispersion: 

BF now < (2/fr)(120/f)(m fl /lGeV)(W~ 3 /v). (4.4) 

Taking v~5x 10~ 4 , and writing = (1/2) (1 + \/5)mo (so that wiq ~ m^, but with some 
additional dependence on 5 and m^), we obtain, 

BF now < (250//)(mo/lGeV). (4.5) 

For example, taking / ~ 0.6 — 0.7, as appropriate for the pure-electron case, yields a maximal 
present-day boost factor of 70-80, for uiq ~ 200 MeV. 

4.1 Other limits from the early universe 

Limits on DM annihilation from the total optical depth to the last scattering surface were 
studied by [52, 53] 5 . These constraints are extremely similar in spirit to the CMB limits 
studied above, measuring the perturbation to the optical depth from the modified ionization 
history, and would seem to be subsumed in the complete likelihood analysis taking into 
account the modified ionization history (which gives rise to the limits quoted above); the 
comparison of the integrated optical depth (due to a nonstandard ionization history) directly 
to the WMAP limit on r (which is derived assuming the standard ionization history) may 



The possibility for DM to contribute significantly to reionization had been previously discussed in [54]. 
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introduce some error relative to the complete likelihood analysis. The limits in [53] appear 
stronger by a factor of ~ 2: however, their approach tends to overestimate (by an 0(2) 
factor) the fraction of deposited energy relative to the careful numerical calculation [38]. 
Limits from the extra heating induced by DM annihilation in the early universe were also 
studied by [53]: they were found to be generically weaker than the optical depth bounds. 

Another possible early-universe limit can be obtained from the extragalactic gamma-ray 
background [52, 55-58], but this constraint tends to depend very strongly on the assumed DM 
structure formation history, and the density profiles of the early halos. Depending on the halo 
density profile and the extrapolation of the halo mass function to low masses, this bound can 
be an order of magnitude either stronger or weaker than the limits derived from the CMB; 
thus, it does not provide any additional robust constraints on the scenarios we consider. 
Furthermore, the DM self-interactions present in models of this type, and velocity "kicks" 
from decays or de-excitations of the excited state, both have the potential to modify structure 
formation in ways which are not accounted for in collisionless cold dark matter simulations 
(see e.g. [59, 60] and references therein). This issue also affects potential constraints from 
DM substructure in the Milky Way. 

4.2 Limits from the Galactic halo 

Measurements of gamma-rays and neutrinos from the Galactic halo have also been used 
to constrain DM annihilation. The gamma-ray signal has two components: (1) photons 
produced directly in DM annihilation by final state radiation, or by production of neutral 
pions which decay into gammas, and (2) the signal produced by inverse Compton scattering 
of high-energy electrons on starlight, infrared and CMB photons. Signals from neutrinos and 
the first gamma-ray component (which we will denote "FSR") depend on the dark matter 
density profile and (because of Sommerfeld enhancement) the dark matter velocity profile; 
the second gamma-ray component (denoted "ICS") has additional astrophysical uncertainties 
associated with the electron propagation. 

The neutrino signal, while somewhat model-dependent (requiring a significant branch- 
ing ratio into non-e + e~ final states), has the advantage of having no significant astrophysical 
background; future studies at IceCube may place very robust constraints on the scenarios 
we consider [61, 62]. However, for e.g. a pure A(i final state, the present limits from Su- 
perKamiokande can only rule out higher-mass DM {m x > 3 TeV) models fitting the cosmic- 
ray excesses, assuming an Einasto DM density profile and no significant change in the DM 
velocity distribution with Galactocentric radius [31]; even accepting these (far from conser- 
vative) assumptions, most of the parameter space we will be interested in lies at lower DM 
masses. 

The gamma-ray bounds [31, 58, 63-69] are often more constraining 6 . However, the 
most recent conservative analyses [58, 68, 69] indicate that DM annihilation models fitting 
the cosmic-ray data remain allowed if the final state consists of electrons and muons, and 
the DM density profile is rather shallow; these models are generally not ruled out by any 
Galactic gamma-ray limits (models with a mixture of electrons, muons and charged pions 
in the final state were not tested in these papers, but since the spectrum is fairly similar 
we do not expect the constraints to change greatly). Shallow DM profiles, like the "cored 

6 It is worth noting, however, that bounds from the H.E.S.S measurement of gamma-ray emission from the 
Galactic Ridge rely on a background subtraction which can also include DM-related emission, and this has 
not always been taken into account in setting limits [70]. 
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isothermal" or Burkert benchmarks commonly used in the literature, are disfavored by DM- 
only N-body simulations, but appear to agree better with observations (see e.g. [71] for a 
recent review) . It is not clear how the inclusion of baryons affects the DM density profile in 
the inner galaxy: it has been suggested that baryons could either steepen the DM cusp via 
adiabatic contraction [72, 73] or flatten it via dynamical friction [74, 75]. 

In Sommerfeld-enhanced models, the uncertainty in the density profile is compounded 
by uncertainty in the variation of the velocity distribution with radius. If the DM velocity 
drops markedly in the Galactic center, as generically expected from DM-only simulations 
(e.g. [76]), the constraints can become more stringent, as the Sommerfeld enhancement 
becomes larger (if it is not yet saturated) ; the shift in gamma-ray constraints in this scenario 
has been considered by [77]. Conversely, if the velocity dispersion of DM particles rises 
in the inner halo, the annihilation cross section will decrease, relaxing the constraints [78]. 
Some simulations including baryons find a rise in the velocity dispersion in the inner galaxy, 
accompanied by a flattened DM core [79-84]. 

Increased tidal disruption is expected to lower the amount of substructure closer to 
the Galactic center (e.g. [85-87]), so if there is a substantial contribution to the PAMELA 
and Fermi signals from DM substructure in the neighborhood of the solar system, this 
could alleviate tension between scenarios fitting the cosmic-ray data and both the CMB 
constraints and the gamma-ray limits from the inner Galaxy [88-91]. Our assumption of 
a Maxwell-Boltzmann distribution for the local DM velocity is also only an approximation 
to the reality: instead employing a velocity distribution based on N-body simulations can 
increase the Sommerfeld-enhanced annihilation rate by up to a factor of ~ 2 [92], which 
improves the range of parameters consistent with the CMB that can also explain the cosmic- 
ray data, and may (depending on the variation in the shape of the velocity distribution with 
Galactocentric radius) also modify the Galactic gamma-ray limits. 

Furthermore, in the scenarios we consider - with their large branching ratios to leptons 
and small direct photon production - the strongest limits tend to come from inverse Compton 
scattering, and astrophysical uncertainties on the cosmic ray propagation are relevant. Most 
of the analyses to date have assumed that essentially all of the electron energy is lost to 
ICS, rather than to synchrotron radiation; however, this assumption is strongly dependent 
on the magnetic field in the inner Galaxy, which is not well constrained (for a recent review 
of B-fields in the Galactic Center see [93] , and [94] for a discussion of fields off the Galactic 
plane). Furthermore, Fermi data indicate the presence of a large-scale structure in gamma- 
rays toward the Galactic center, extending roughly 50° north and south of the plane [95]. 
It has been suggested that this structure could be associated with an AGN jet or Galactic 
wind; if this is the case, cosmic ray propagation in this region - which includes the region 
used to set some of the most stringent constraints - could be greatly modified, and ICS limits 
from this region assuming steady-state diffusive propagation could be very inaccurate (for 
example, electrons - both background cosmic rays and DM annihilation products - may be 
swept away by the jet or wind before they can scatter). 

4.3 Limits on the force carrier mass 

Measurements of gamma-ray emission from dwarf galaxies have been used to constrain 
Sommerfeld-enhanced models specifically, since the velocity dispersion in these structures 
is thought to be lower than in the local halo [96]; at present, mediator masses < 100 
MeV are disfavored by these constraints. Constraints on the dark matter self-interaction 
cross section have been discussed recently by [35, 97], with the result that mediator masses 
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m <t> ^30 — 40 MeV can be excluded. Consequently, we choose all our benchmark mediator 
masses to exceed 100 MeV. 

5 The Parameter Space for PAMELA and Fermi 

With the results of the previous section, we can now answer the question: for what, if any, 
parameter choices can the cosmic ray signals of PAMELA and Fermi be explained, while 
still achieving the appropriate relic abundance and evading the CMB limits? Determining 
whether such regions of parameter space exist can be somewhat involved, because of the 
interplay between the various parameters. Specifically, the particular Sommerfeld enhance- 
ment depends greatly on m^, but so too does the decay mode of (f>, and thus the needed cross 
section as well. As such, we approach this in a systematic fashion. First, for a representative 
set of decay modes for eft, spanning a range of masses = 200 — 900 MeV, we determine the 
appropriate ranges of cross section and mass that can explain the data. Having found the 
preferred values, we then study the Sommerfeld enhancement along slices of fixed and 
m x , keeping the relic abundance fixed to the observed value. Finally, we select "benchmark 
points" in this parameter space to give representative examples of explicit fits to the data. 

To determine whether a particular choice of (m x ,m^,'BF) explains the e + e~ signals, we 
consider a broad range of parameters describing the background e + e~ spectra and cosmic ray 
propagation. We use the publicly available GALPR0P code [98] to calculate the local cosmic ray 
(CR) signals as measured by Fermi and PAMELA. GALPR0P allows for considerable freedom 
in the details of cosmic ray production and propagation. Among other things, one may choose 
the electron and proton injection spectra, i.e. the spectra of primary electrons and protons 
emitted at the source; the diffusion parameters, which include the diffusion coefficient and 
the size of the cylindrical diffusion zone; and the Galactic magnetic field and interstellar 
radiation field, which determine cosmic ray energy losses. 

5.1 Background cosmic ray spectra 

The background of CR electrons consists of "primary" and "secondary" electrons; primary 
electrons are injected into the Galaxy by cosmic ray sources such as supernova remnants, while 
secondary electrons are created in the interactions of CR protons with the interstellar gas. 
The background of positrons is entirely secondary in nature, arising from the interactions of 
the CR protons with the interstellar medium. Therefore, both the primary electron spectrum 
and the primary proton spectrum are fundamental components of our analysis of the fits to 
PAMELA and Fermi. Both the proton and electron injection spectra are broken power laws in 
energy and can be described by four parameters each (n, 71, 72, Eb), an overall normalization 
n, the power law indices below and above the break 71 and 72, and the energy at which the 
break occurs E^ . 

7 The GALPR0P code describes the primary electron spectrum by a double broken power law in energy, 
which requires two additional parameters, a second break energy and a third power law index. The spectrum 
typically has two breaks above ~ 100 MeV, one around a few GeV and a second around 1 TeV. The low 
energy break results from the energy dependence of the dominant energy loss processes changing from oc 1 /E 
to oc 1/E 2 . A high energy break is expected to exist due to the finite distance CR electrons can travel from 
the source before losing a large portion of their energy. The distances to the nearest CR sources, for example 
supernova events, determine the highest energy the background electrons reaching our detectors can have. 
Since our PAMELA and Fermi fits do not extend to energies below 10 GeV, the low energy electron spectrum 
is irrelevant. Therefore, for our purposes, the primary electron spectrum can be described by a broken power 
law with two indices and one break energy. 
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We constrain the primary proton spectrum by requiring the propagated proton spectrum 
to be in agreement with the local proton measurements by AMS-01 [2], BESS [99], and IMAX 
[100] above 1 GeV, and the propagated antiproton spectrum and the antiproton-over-proton 
ratio p/p to be in agreement with the recent PAMELA measurements [101, 102]. Accordingly, 
for the protons we take the power law index below E\, = 9 GeV to be 71 = —1.98. We note 
that our fits to the PAMELA and Fermi e + e~ data do not extend below 10 GeV, so the 
spectrum of protons below 10 GeV does not affect our fits. Therefore, we do not vary the 
form of the proton spectrum below the break energy. Above 9 GeV, we allow the power law 
index to vary such that the resulting spectrum does not exceed the local data by more than 
3a. The PAMELA p/p data place the strongest constraint on the range of acceptable power 
law indices for the injection spectrum of primary protons above 9 GeV. We find the allowed 
range to be -2.50 < 72 < -2.10. 

We use the Fermi [6] electron data to constrain the injection spectrum of primary 
electrons above 10 GeV and confirm our findings with several other measurements of the 
electron spectrum [103-106]. Fitting to the Fermi data between ~ 70 GeV and ~ 380 GeV, 
we find that the range of indices —2.57 to —2.18 gives spectra that are consistent with the data 
such that no Fermi data point is exceeded by more than 3a. We stress that these constraints 
arise from fitting only the primary electrons to the e + e~ Fermi data; we have not included 
the fluxes of secondary electrons and positrons in our fits. Since the secondaries are an order 
of magnitude or more smaller than the flux of primary e~'s between 10 GeV and 100 GeV, 
including them would result in only a small softening of the allowed spectrum, making the 
range of indices slightly more negative. However, inclusion of a dark matter component of 
e + e~ in the flux allows for considerable softening of the primary electron spectrum. To 
accommodate the effects of a large dark matter e + e~ component, we extend the range of 
allowed indices to smaller values and perform our scans over the range —3.00 to —2.00. Our 
fits to the PAMELA and Fermi data confirm that this is a reasonable choice, as only a few 
of the 2a Fermi fits require an electron injection spectrum with an index smaller than —3.00 
and none of the fits require an index larger than —2.40. The H.E.S.S data indicate a steep 
falloff in the e + e~ spectrum above ~ 1 TeV [7, 107]. We assume that the local primary 
electron spectrum exhibits this behavior and introduce a break at 2.2 TeV, above which the 
injection spectrum has a power law index of —3.30. 

5.2 Propagation parameters 

We define the cylindrical diffusion zone centered on the Galactic Plane (GP) to have a radius 
of = 20 kpc and height hd = ±4 kpc. While one may find reference to diffusion zone heights 
as small as hd = 1 kpc and as large as hd = 15 kpc, these are extreme values. Standard values 
lie in the range 3 — 7 kpc. Inside the diffusion zone, cosmic ray propagation is governed by the 
diffusion-loss equation. At the boundaries of the diffusion zone, free escape of all cosmic rays 
is assumed to occur. Measurements of the local value of the unstable secondary-to-primary 
ratio Be 10 /Be 9 provide a good check on the choice of diffusion zone height, since this ratio of 
secondaries-to-primaries is dependent on the average time cosmic rays spend in the Galaxy 
before free-escape, decay, and spallation. The energy dependent diffusion coefficient has the 
form D(E) = D x 10 28 (E/A GeV)° cm 2 s _1 . Typical values of a are between 0.33 and 
0.50 (with a = 0.33 corresponding to a turbulent magnetic field with a Kolmogorov type 
spectrum and a = 0.50 corresponding to a turbulent field with a Kraichnan type spectrum), 
while Dq takes values in the range ~ 2 — 10. The stable secondary-to-primary ratios B/C and 
sub-Fe/Fe are a measure of the average amount of matter cosmic rays propagate through. For 
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this reason, measurements of B/C and sub-Fe/Fe help constrain the diffusion coefficient. We 
calculate our fits to the PAMELA and Fermi data using the benchmark diffusion parameters 
-Do = 4.00 and a = 0.50, which are broadly similar to the Ml parameters of [108], albeit 
with a more conventional disk scale height (i.e. hd = 4kpc rather than the 15 kpc of the 
Ml model), and fit the data for these three secondary-to-primary ratios at the 3a level. The 
effect on the local cosmic-ray spectra of variations in the diffusion parameters can largely 
be compensated for by changes in the background spectral indices, so we choose to fix the 
diffusion parameters and vary the background as described above. 
We model the magnetic field as an exponential disk, 



(fl a > 1 / a = Bbexp 



(r-R&) \z\ 

TB Z B 



(5.1) 



where r# = 10 kpc is the radial scale, zb = 2 kpc is the vertical scale, and Bq = 5.0 (iG is 
the local value. Additionally, we use the most current model of the interstellar radiation field 
available to the public, that of [109]. 

5.3 Dark matter density 

For our choice of dark matter density profile, we use an Einasto profile [110, 111] 

2 fr a -R° r . 



p(r) = po exp 



a 



(5.2) 



with a = 0.17 and r_2 = 25 kpc. The local density is uncertain, but the best current estimates 
are p /(GeV cm" 3 ) = 0.385 ± 0.027 [112], 0.43 ± 0.11 ± 0.10 [113], and 0.466 ± 0.033 ± 0.077 
[114]. We conservatively take po = 0.4 GeV/cm 3 , but note that taking a ~ la high value 
from the result of [114] would lower the needed boost by a factor of 2. Since the majority of 
the e + e~ signal is produced in the nearby ~ 1 kpc, the question of the dark matter density 
profile - Einasto, NFW, cored isothermal or something else - is not particularly important. 

5.4 Quality of fits 

Others have previously employed x 2 analyses to determine the best-fit regions for particular 
final states [10, 31, 115]; however, in the absence of a covariance matrix for the data, the 
information that can be extracted from the x 2 is limited. The usual assumption is that 
the error bars are completely uncorrelated, and this is indeed a conservative approach for 
the simple question of whether a given model is a good fit (x 2 /dof < 1 + e). However, 
the question of whether one model with \ 2 1 is a better fit than another, and with 
what confidence, cannot be answered without the covariance matrix, since if one simply 
assumes uncorrelated errors then the data is skewed by repeatedly counting correlated errors, 
erroneously disfavoring regions of parameter space in an unpredictable way. 

As an example, imagine we have a dataset with perfectly correlated large systematic 
errors (~ a), and very small statistical errors (~ e). Suppose we have two models, one of 
which is consistently ~ la above the data, and one of which varies smoothly from ~ la 
above the data at the lowest energies to ~ la below the data at the highest energies. If we 
incorrectly assume that all the errors are uncorrelated, and add the systematic and statistical 
errors in quadrature, then the x 2 /dof will be ~ 1 for the first model and ~ 1/3 for the second 
model; if the number of degrees of freedom is large, a naive Ax 2 analysis will find that 
the second model is preferred at very high confidence. If we take into account the perfect 
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correlation of the systematic errors, however, then the first model still has x 2 /dof ~ 1, but 
the x 2 /dof for the second model balloons by a factor of ~ (cr/e) 2 . The conservative approach 
of taking all the errors to be uncorrelated will never rule out a model that is actually a good 
fit, but a model that appears to be the best fit may no longer be the best fit, and in fact may 
be ruled out, once the full covariance matrix is taken into account. 

Consequently, to determine if a particular choice of (m x ,m^,6) fits the data, we insist 
that the background+signal curves pass through the error bars for the PAMELA positron 
fraction data above 10 GeV and the Fermi e + + e~ data above 30 GeV (where charge- 
dependent and charge-independent solar modulation, respectively, are thought to be small), 
for some set of parameters governing the background electron and positron spectra, within 
the bounds described in §5.1. We compute the regions where these conditions are satisfied 
for the la, 90% confidence, and 2a error bars. We additionally require that the curves not 
exceed the H.E.S.S e + + e~ systematic error band, which is, in effect, forcing the e + + e~ 
spectrum to satisfy the E~ A fall-off above ~ 1 TeV. Using the la errors, these requirements 
are more constraining than demanding x 2 /dof < 1 + e, allowing us to select a preferred 
region in the allowed parameter space without having to explicitly write down a model for 
the covariance matrix. This approach has the positive features that: 

• It disfavors models where a small number of points lie well outside the error bars; since 
the true uncorrelated fluctuations are likely to be much smaller than the nominal error 
bars, this is a positive feature. 

• It treats different datasets on an equal footing; the favored parameters provide model 
curves in good agreement with PAMELA and Fermi. In contrast, with a x 2 test 
assuming uncorrelated errors the "best-fit" models tend to have x 2 /dof <C 1 for the 
Fermi data, which compensates for a somewhat poor fit to the PAMELA data. This 
skews the "best-fit" regions towards lower values of the boost factor. 

This method makes no claims to determine the relative likelihoods of models that each fit the 
data well. In the absence of covariance matrices from the PAMELA and Fermi experiments, 
it is probably the most reasonable quantitative test of "explaining the data" that is available 
to us. 

With regard to the H.E.S.S data, two points are relevant here. First, we follow [116] 
and treat the H.E.S.S data as an upper limit. Although photon contamination is expected 
to be low at these energies, the data are still properly limits. Second, as noted by [116], the 
data themselves seem inconsistent with the Fermi data. Attempting to fit both is then likely 
to be partially fitting to a systematic uncertainty. Consequently, we do not fit the H.E.S.S 
data, but do show them for our benchmark points, rescaled to be consistent with Fermi. We 
differ slightly from [116] here, in that we assume that a small (~ 8%) shift in the energy scale 
explains this, rather than a 15% shift in overall normalization. While these data are not 
explicitly included in our fits, our benchmark points nonetheless agree well with the rescaled 
H.E.S.S data. 

As a check on our method, we also computed the x 2 as a function of ?tidm and BF, 
choosing the best-fit background parameters at each point. The allowed regions - that is, the 
regions that are not excluded by a goodness-of-fit x 2 test, at 90%, 95% and 99% confidence 
— are quite large, and as mentioned above, they tend to skew toward lower boost factors due 
to the large number of Fermi data points (with covariances that are not taken into account). 
Inclusion of the rescaled H.E.S.S data in the \ 2 (instead of using H.E.S.S only to set upper 
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limits) substantially reduces the size of the allowed regions; however, all of our benchmark 
points, and the bulk of our preferred parameter space, remain allowed. Some small regions 
with high boost factors and relatively low DM masses are ostensibly ruled out: this is because 
such models require a fairly soft background electron spectrum, to fit the data with a large 
DM contribution, and at energies above the mass of the DM this soft background spectrum 
undershoots the H.E.S.S data. It is possible that by varying the position of the high-energy 
break in the background electron spectrum, or allowing a small shift in the overall energy 
scale, a better fit might be achievable even for these models, but in any case, the effect of a 
conservative \ 2 analysis is generally to shift the preferred range of boost factors downward. 

6 Results 

6.1 Parameter scans 

To determine the current annihilation rate, we scan a three-dimensional parameter space 
consisting of the DM mass m x , the dark gauge boson mass m<*, and the mass splitting 5. 
We focus on the singly-charged Higgs model with the mass splitting generated as described 
in §2.2, which we refer to as the "minimal model". The dark sector coupling ao is deter- 
mined by requiring the thermal relic density to be O/i 2 = 0.1120, in accordance with the 
WMAP7+BAO+H0 ML cosmological parameter set [117-119]. The fits to the PAMELA 
and Fermi cosmic-ray excesses are determined by m x , which sets the overall energy scale, 
and m^, which determines the branching ratios of the dark gauge boson to different SM final 
states (as summarized in e.g. [120]). 

We show in Figure 1 the allowed BF-mass parameter space for four representative final 
states, corresponding to four benchmark mediator masses: cp — > e + e~ (m^ = 200 MeV), 
(f> -> e+e" :<j) ->■ /x+/x~=l:l (m<p = 350 MeV), (f> -> e+e" :</) ->■ fi + ^:<f> -> 7r+7r~=l:l:l 
(m^ = 580 MeV), and <f> -> e+e~ :</> -> ^ + ^:<j) ->■ 7r+7r~=l:l:2 (m^ = 900 MeV). At higher 
mediator masses, hadronic final states become increasingly important and the softness of the 
spectrum makes it more challenging to fit the data. In the left-hand column of Figure 1 we 
show the fits lying within the la, 90% confidence, and 2a error bars for PAMELA only, for 
Fermi only, and for both PAMELA and Fermi. In the right-hand column we again show the 
fits within the la, 90% confidence, and 2a error bars for both PAMELA and Fermi. 

We briefly clarify some of the assumptions made in determining the parameter space 
shown in Figure 1. As discussed in §5.1 and §5.2, we describe the background proton spectrum 
as a power law with one break. For our analysis, we fix both the break energy and the index 
below this break energy. See Table 1 for values. We allow the upper index 7 P 2 to vary in 
the range (—2.50, —2.10) and the overall proton normalization to vary. Similarly, we describe 
the electron spectrum as a power law with two breaks. We hold the spectral index fixed 
for energies below E e \ = 4 GeV. Moreover, the high-energy index, 7 e 3 = —3.30, and break 
energy, E e 2 = 2200 GeV, are held fixed at values chosen to give agreement with the Fermi 
and H.E.S.S data at energies above ~ 1 TeV. We are left with two varying parameters for 
the electron spectrum, the electron spectral index in the energy range 4 to 2200 GeV, which 
is allowed to vary between —3.00 and —2.00, and the overall normalization of the spectrum. 
In summary, the BF-mass parameter space shown in Figure 1 is determined by varying four 
parameters, the proton spectral index above 9 GeV (which can take values in the range 
(—2.50,-2.10)), the electron spectral index between 4-2200 GeV (which can take values in 
the range (—3.00, —2.00)), and the normalization of each spectrum. As discussed in §5.2, the 
effects of the diffusion parameters on the background cosmic-ray spectra can be mimicked 
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by a change in the spectral indices. Thus, a variety of diffusion scenarios are represented in 
the analysis leading to the results shown in Figure 1. A BF-mass combination is considered 
to fit the data if, within these constraints, the background parameters can be chosen so 
the signal+background curve passes through every error bar (for example, the "lcr" regions 
require the model curve to pass through every la error bar). 

We include curves in Figure 1 describing the allowed values of the local boost factor 
(BF) for several different DM mass splittings 5, subject to the relic density constraint. Recall 
that the boost factor is defined as (av}/3 x 10 _26 cm 3 s _1 . The solid portions of the curves 
are those for which the CMB constraints are met. In both sets of plots we mark with a "+" 
benchmark points for each mediator mass; a benchmark point is one that simultaneously 
fits the PAMELA and Fermi data with the correct relic density, while satisfying the CMB 
constraints. Results are shown for 800 GeV < m x < 3 TeV only, though the allowed regions 
may extend to lower and higher masses. For comparison to earlier work (e.g. [33, 34]), we 
also show the case with no mass splitting; in this case we omit annihilation channels involving 
the dark Higgs, since such annihilation channels were not included in the previous studies. 
It is clear that the introduction of an O(MeV) mass splitting removes any tension between 
the preferred PAMELA/Fermi region and the boost factor from Sommerfeld enhancement. 
However, the CMB constraints introduce significant tension at lower mediator masses. 

In Figure 2, we hold the DM mass fixed and scan over and 5, for m x = 0.9, 1.2, .1.5 
and 1.8 TeV in the minimal model. For illustration, for m x = 1.2 TeV we also show the 
effect of changing the coannihilation parameter to k = 1,4. For k > 1/4 we include only the 
self-annihilation channel XX ~~ ^ <^> and a coannihilation channel parameterized by k; with 
these assumptions the only dependence of the annihilation rate on and 5 comes through 
the low-velocity Sommerfeld enhancement, so except near the resonance centers, the relic 
density is largely fixed by ao and m x . Since these resonance regions are generally already 
ruled out by constraints from the CMB, we simply hold olb fixed for these scans, and confirm 
that an acceptable relic density is obtained in all the regions that are not ruled out. For the 
minimal model, in contrast, the early- universe cross section for the new annihilation channel 
XiXi hoho scales as <5 2 /m^, so the coupling an must be reduced at small and/or large 
8 to obtain the correct relic density. 

A quick study of these plots shows that boosts larger than 200 populate a large region of 
the allowed parameter space, with boosts even larger than 300. Much larger boosts are still 
possible to achieve while maintaining agreement with the relic density constraint (S ~ 500) 
but are strongly disfavored from CMB constraints (note also that not every point with a high 
BF will provide good fits to the CR data; in particular, the p resonance leads to a pronounced 
dip in the leptonic fraction for ~ 0.8 GeV that can make the spectrum too soft to provide 
a good fit). For low or large 5, the XX ~~ ^ hh annihilation channel in the minimal model 
naturally reduces the annihilation rate both in the local halo and during the epoch of last 
scattering; increasing k has the same effect, but independent of and 5. Either opens up 
new regions of parameter space at low mediator masses, which would be ruled out by the 
CMB constraints in the minimal model with 5 = 0. 

6.2 Benchmark points 

We now provide specific examples of benchmark points that satisfy the relic density and 
CMB constraints. Table 1 lists the particle physics parameters, present-day and low-velocity 
boost factors, and the CMB limit on the boost factor for these points. The boost factors 
are computed using the semi- analytic approximate formula for the Sommerfeld enhancement, 
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Figure 1. Left: Allowed ranges of parameter space for fits within the la, 90% confidence, and 2a error 
bars to PAMELA only (in decreasing intensity of red), Fermi only (in decreasing intensity of gray), and for 
simultaneous fits to both PAMELA and Fermi (in decreasing intensity of purple). Yellow crosses indicate 
benchmark points. Right: As in left, with curves showing the boost factors for a range of mass splittings 8 such 
that Qh 2 = 0.1120 (dashed). Yellow lines, marked with asterisks, are chosen to pass through the benchmark 
points for cases where the BF varies rapidly with S. The CMB constraints are met for the solid portions of 
the curves. Results are shown for 800 GeV < m x < 3 TeV only. All preferred regions shown here assume 
po = 0.4 GeV/cm 3 and no contribution to the signal from DM substructure; any substructure correction (e.g. 
[87]) will shift the preferred regions to lower boost factors. The (5 = curve is intended as a consistency check 
with previous work, and so annihilation channels involving the dark Higgs were omitted from the computation 
in this case. 
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m (GeV) m (GeV) m, (GeV) 

Figure 2. Contours for the boost factor BF in the local halo as a function of the mediator mass 
and mass splitting 8, with an chosen to produce the correct relic density. The DM velocity distribution is 
assumed to be Maxwellian with a = 150 km/s. The regions overlaid with red lines are ruled out at 95% 
confidence by bounds from WMAP5, taking / = 0.2 for hadronic final states, / = 0.24 for muons, and / = 0.7 
for electrons, and weighting these three contributions according to the (m^-dependent) branching ratios for 
decays of the dark gauge boson. Blue-hatched regions illustrate the effect of shifting the WMAP5 constraint: 
these parts of parameter space would remain ruled out even if the bounds were relaxed by a factor of 4 (e.g. 
due to degeneracy with some parameter not included in the analysis, although we do not consider this a likely 
scenario). Upper row: Results for 1.2 TeV dark matter, for three values of the coannihilation parameter k. 
Upper left panel: k = 1/4, corresponding to the minimal singly-charged Higgs model described in §2.2. Upper 
center panel: k = 1 (ao = 0.0263, see text). Upper right panel: k — 4 (old = 0.0177, see text). Lower row: 
Results for DM mass (lower left panel) 900 GeV, (lower center panel) 1.5 TeV, and (lower right panel) 1.8 
TeV, in the k = 1/4 minimal model. Capture into a bound state, inducing an additional enhancement to 
late-time annihilation, is kinematically allowed in regions to the left of and below the black-dashed line, but 
has not been included in the analysis; see Appendix A for a discussion. 
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Benchmark 


Annihilation 




m x 




r 



Local 


saturated 


OlvlhS 


number 


channel 


(MeV) 


(TeV) 




(MeV) 


BF 


BF 


limit 


1 


1:1:2 e* : : tt* 


900 


1.68 


0.04067 


0.15 


300 


530 


600 


2 


1:1:2 : (j* : tt* 


900 


1.52 


0.03725 


1.34 


260 


360 


545 


3 


1:1:1 e ± :fj, ± : tt* 


580 


1.55 


0.03523 


1.49 


250 


437 


490 


4 


1:1:1 : (j* : tt* 


580 


1.20 


0.03054 


1.00 


244 


374 


379 


5 


1:1 : ^ 


350 


1.33 


0.02643 


1.10 


156 


339 


340 


6 


e only 


200 


1.00 


0.01622 


0.70 


67 


171 


171 



Table 1. Particle physics parameters, present day boost factors, and boost factors in the low- velocity 
limit for benchmark points. The boost factor (BF) is defined as (av)/3 x 10 -26 cm 3 s _1 . 



Benchmark 


n e X 10~ iu 


7el 


E e i 


7e2 


E e 2 


7e3 


n p x 10 a 


7pi 


E p \ 


7p2 


number 


@ 34.5 GeV 




(GeV) 




(GeV) 




@ 100 GeV 




(GeV) 




1 


3.12763 


-1.60 


4.0 


-2.45 


2200 


-3.3 


5.66361 


-1.98 


9.0 


-2.11 


2 


3.12763 


-1.60 


4.0 


-2.45 


2200 


-3.3 


5.66361 


-1.98 


9.0 


-2.11 


3 


3.12763 


-1.60 


4.0 


-2.45 


2200 


-3.3 


5.66361 


-1.98 


9.0 


-2.11 


4 


3.06444 


-1.60 


4.0 


-2.50 


2200 


-3.3 


5.20440 


-1.98 


9.0 


-2.11 


5 


3.09604 


-1.60 


4.0 


-2.45 


2200 


-3.3 


5.74014 


-1.98 


9.0 


-2.11 


6 


3.10299 


-1.60 


4.0 


-2.45 


2200 


-3.3 


5.74014 


-1.98 


9.0 


-2.11 



Table 2. GALPROP parameters describing the electron and proton injection spectra for the benchmark 
points. For all cases the diffusion parameters used are the following: Dq = 4.00 (multiplied by 
10 28 cm 2 s _1 to obtain D(E) at E = 4 GeV), a = 0.50, hd = 4.0 kpc. Normalizations n e and n p are 
in units of cm -2 s _1 sr _1 MeV" 1 . Because we do not fit the PAMELA or Fermi data below 10 GeV, 
the values of E e \ and 7 e i are unnecessary for our fits. However, we list our default values here. 



but numerical checks indicate that the approximation is accurate to within ~ 5% for v > 150 
km/s, and to within ~ 10% for v — > 0. In both cases, the approximation tends to overestimate 
the enhancement, so the overall effect is to slightly weaken the CMB constraints relative to 
the present-day boost factor. 

Two of our benchmarks have mass splittings well below 2m e , and therefore potentially 
long lifetimes. We can estimate the depletion of the excited state due to DM-DM downscat- 
tering using the prescription for the scattering rate described in §3.1, although it should be 
noted that unlike the relic density calculation, the uncertainties in our prescription for the 
downscattering cross section induce correspondingly large uncertainties in the relic excited 
fraction. For the 1.68 TeV and 1 TeV benchmarks, the estimated relic excitation fractions (af- 
ter decoupling of the DM-DM downscattering, but before decay) are respectively ~ 5 x 10~ 3 
and 2 x 10" 3 . 

In Table 2 we provide the GALPROP parameters for the electron and proton injection 
spectra needed to reproduce the background e + e~ spectra for each of the benchmark points. 
As discussed in §5.1, the injection spectrum for electrons is a power law in energy with 
two breaks and so can be described by six parameters (n e , 7 e i, E e i, j e 2, E e 2, 7e3), an overall 
normalization n e , the three power law indices, j e \, j e 2, and 7 e 3, and the energies at which 
the breaks occur E e \ and E e %. The low energy break occurs around a few GeV. Since our 
PAMELA and Fermi fits do not extend to energies below 10 GeV, the electron spectrum at 
these energies is irrelevant for our analysis. Therefore, the parameters E e \ and 7 e i are un- 
necessary for our fits, but we include them in Table 2 for completeness. The proton injection 
spectrum is a broken power law in energy described by the four parameters (n p , ^ p i,E p , 7^2)- 
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Figure 3. Benchmark models fitting the PAMELA (first and third rows) and Fermi (second and fourth 
rows) cosmic-ray excesses, obtained using the GALPROP program. 



6.3 Comparisons with Previous Results 

While we find ample regions of parameter space that provide agreement with the PAMELA 
and Fermi results, previous studies [35] have been more negative. In particular, [35] finds 
a maximum local "boost factor" (BF) of ~ 120 from Sommerfeld-enhanced annihilation for 
~ 2 TeV DM, and a BF of 90 for 1 TeV DM, compared to a best fit to the data for XX ~~ > 
with 4> -> n+yT (taken from [31, 116]) of 2.35 TeV DM with a boost factor of 1500. We 
should emphasize that our results are completely consistent with theirs, with the different 
conclusions arising from our consideration of a more general parameter space. Specifically, 
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• A decay mode of <p — > p + p~ was assumed by [35], which is natural for models where 4> 
is a scalar, but does not occur in models where eft is a vector, and the force arises from 
a conventional gauge group. In these cases, (p couples to charge, and there is always 
a sizeable hard <j) —> e + e~ component, unless <j) is very degenerate with the p meson. 
The presence of an electron component hardens the final e + e~ spectrum, increases the 
power in e + e~ as opposed to neutrinos (by a factor of up to ~ 3, depending on the 
branching ratio), and lowers the preferred mass scale for the DM to around 1 TeV, 
since the electron component dominates the high-energy cutoff behavior which sets the 
mass scale, both of which lower the needed boost. 

• For vector mediators, it is both natural and almost an experimental necessity (given 
constraints from direct detection) to consider the case where xi an d Xi are non- 
degenerate. Such a scenario generally produces larger Sommerfeld enhancements than 
in the degenerate case. 

• Older works on the local DM density [121] found a central value of po = 0.3 GeV/cm 3 , 
a number which has become a standard. As we noted above, more recent studies, 
however, find values of 0.39 [112], 0.43 [113] and 0.46 GeV/cm 3 [114]. Taking the 
current best estimates of the local density then leads to a factor of ~ 2 reduction in 
the required boost factor. 

The combination of these effects suggests that the natural boost required for the major- 
ity of the parameter space is 0(100-300), or a factor of 0(5) lower than the models considered 
in [35], or an order of magnitude (relative to [31, 116]) when combined with the best current 
estimates for the local relic density. Our results are consistent with the conclusions of [35] 
that XX </></>> 4> ~ > is tightly constrained, but the general parameter space (in which 

some <j) —> e + e~ is present) is far less so. This conclusion need not invoke artificially high 
local substructure, but only a more accurate estimate of the preferred boost factor and DM 
mass range for the specific class of scenarios we consider, where the DM annihilates through 
a dark gauge boson that kinetically mixes with SM hypercharge. 

7 Conclusions 

Data from several experiments have pointed to the presence of a new, primary source of 
cosmic ray electrons and positrons. Dark matter is a long-standing candidate for such a 
source, but has difficulty achieving the high rates, hard spectrum and dearth of correlated 
anti-protons that are required by the data. Models of dark matter with a new, light boson 
(p can qualitatively explain such phenomena through annihilations XX ~ > </></>) followed by 
<j) — > e + e~ , p + p~ , ir + ir~ . 

We have seen that this agreement extends to fully quantitative connections as well. 
Significant regions of parameter space exist for general models where the present-day boost 
is large enough to explain the observed cosmic-ray signals, while yielding the appropriate 
relic density and evading CMB constraints. We have additionally specified benchmark points 
which are representative of a sizeable fraction of the parameter space, in order to give a first 
set of parameters to be searched for in terrestrial experiments. 

In general, the allowed regions of parameter space have a non-trivial 4> — > e + e~ branch- 
ing ratio - models that go exclusively into p + pT have difficulty in achieving the necessary 
boost while being consistent with other constraints, although a proper treatment of capture 
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into bound states can alleviate this tension for DM masses above 2 TeV. At the same time, 
models that go exclusively to e + e~ are more severely constrained by limits from the CMB, 
as mediators lighter than 2m^ often yield too large a signal in the epoch of recombination. 

This suggests a preferred mass range of 2m „ < ~ lGeV. 

Fortunately, this range of parameters is very testable. Low-energy experiments [122- 
126], such as APEX should be able to study a wide range of parameters, while LHC and 
Tevatron searches [127, 128] offer complementary reach. Finally, for all of these models, the 
Planck satellite should see a signal in its polarization spectrum. In light of this, these models 
remain a viable, testable and thus exciting scenario for dark matter. 

As a resource for interested readers, a web application for the calculations given in this 
work can be found at http://astrometry.fas.harvard.edu/mvogelsb/sommerfeld/. It 
allows the user to compute the present-day and saturated boost factors, the relic density, and 
an estimate for the relic abundance of the excited state (prior to any decay), for our minimal 
model with any choice of parameters (a, m x ,ms, S), with the option to include WIMPonium 
formation in the 5 = case. It also automates the calculation of the semi-analytic approxi- 
mation for the Sommerfeld enhancement, for a U(1)d model, with the initial-state particles 
in either the ground or excited state. Readers wishing to run large parameter scans or ex- 
plore options not available in the web application should contact TS at tslatyerOias . edu 
for access to the underlying code. 
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A Capture into WIMPonium and Heavy DM Scenarios 

The usual Sommerfeld enhancement calculation neglects diagrams describing the capture of 
two DM particles into a bound state - referred to as WIMPonium, in analogy to positronium 
- accompanied by the radiation of a dark gauge boson. This process is kinematically allowed 



which at low velocities, T < m^, reduces to < c? D ra x j\ as stated in e.g. [12]. 

Above the symmetry breaking scale, the <f> is massless and this process is exactly anal- 
ogous to positronium formation with the replacement m e — > m x , «ew ~~ * a D- The non- 
relativistic cross section (valid in the regime < T < m x ) is given approximately by 




(A.l) 
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[129], 

' l£ / u \ {™n\ ( Trap/v \ ( (a D /2v) 2 \ 3 
3 \a 2 D m x /A) \m^U J \l- e -^ D /v ) \l + {a D /2v) 2 J 




v > a D /2, (A.2) 
, v < a D /2, 



where to is the energy of the radiated 0, uj ~ a D m x /A + m x v 2 , again assuming b<1. This 
expression agrees with [12] in the low-velocity limit when is set to zero. The capture cross 
section experiences Sommerfeld enhancement just as the annihilation cross section does, but 
relative to the annihilation cross section is greatly suppressed at high velocities, with av 
scaling as (2v / ao)~ 4 for v/ap > 1/2. At low velocities, on the other hand, the capture cross 
section scales exactly as the Sommerfeld-enhanced direct annihilation. 

This behavior can be understood from the relative ranges of the various interactions: 
the length scale associated with the annihilation operator is 1 /m x (determined by the mass of 
the x)> whereas for the capture operator it is \ jar>m x (determined by the Bohr radius of the 
WIMPonium). The momentum of the incoming particles is m x v: for nonrelativistic particles, 
this is never large enough to probe the r < l/m x region relevant for annihilation, but it is 
large enough that for v > old, capture cannot be treated as a contact interaction. For v < ao, 
on the other hand, both capture and annihilation behave as contact interactions, and so in 
both cases the Sommerfeld enhancement depends only on the behavior of the wavefunction 
at the origin. At least in the elastic (5 = 0) case, the shape of the wavefunction at the origin 
is completely determined by the requirement that it be finite, with the longer-range physics 
only setting its amplitude: consequently, the enhancement from the long-range interaction 
must scale in the same way for any operators localized at the origin. At energies below 
m^, therefore, we employ the expression for WIMPonium capture given by [12], but with 
the replacement irao/v — > S, where S is the Sommerfeld enhancement due to the Yukawa 
potential: 

Wv rcl ) = — — y-^- j I — z- 15, v 4> = ^l- {Am 4> la 1 m x y. (A.3) 

One might ask whether it is justified to neglect WIMPonium capture as we have done 
in previous sections. If > a D m x /4, i.e. the capture process is kinematically forbidden 
at low velocities, then v ~ old/2 corresponds to T ~ a D m x /4 < m^. Thus the strong 
(2v/aD)~ 4 suppression of the capture cross section holds for all temperatures above the 
symmetry breaking scale. 

Furthermore, for T > m^, interactions with the bath of relativistic $>'s can dissociate the 
positronium [12], so the timescale for a WIMPonium bound state to undergo ^-dissociation 
must be compared to the annihilation rate. Only at T < can the capture cross section 
be used as a proxy for the annihilation cross section, since any particles that form a bound 
state will eventually annihilate. 

The lifetime of positronium is r ~ 1.2 x 10~ 10 s for parapositronium (1/4 of the total 
positronium formed) and r k 1.4 x 10~ 7 s for orthopositronium (3/4), scaling as l/m e a| w 
and l/m e ctg W respectively 8 . The lifetimes for the analogous WIMP states are shorter by 
a factor of (m e /m x )(aEW / a D) 5 ' e ■ The photodissociation cross section for positronium is 

8 The lifetime of ortho WIMPonium can be much shorter than expected from this scaling relation - com- 
parable to that of paraWIMPonium - if the <j> can decay into other dark-sector states, due to the availability 
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given by a ~ 2.5 x 10~ cm 2 at threshold, falling roughly as energy -3 above threshold, and 
scaling as oewOo ~ l/ a EW"ie [130]; for an initial estimate, we will just use the threshold 
value, since at earlier times the increased 4> density will cancel out the reduced cross section. 
The 4> density is obtained from Equation 3.2, and can be written n « lO 4O (m0/lGeV) 3 
cm -3 when T ~ mj,. 

Photodissociation and decay are equally likely at T ~ txia when n^av ~ 1/r: this 
relation is equivalent to the conditions, 

m * ~ 1, paraWIMPonium, (A.4) 



a 2 D m x /4 



-1/3 



0.1 , orthoWIMPonium. (A.5) 



2 D m x /4 ~ ' Uo / 



Thus photodissociation is efficient at suppressing paraWIMPonium decay for 3> a^m x /4 
(the ratio of the rates near threshold scales as mV), and always dominates orthoWIMPonium 
decay unless <C a 2 D m x /4. 

The combination of this effect with the suppressed high-velocity cross section justifies 
our neglect of WIMPonium capture for > a 2 D m x /A, which holds true for all our bench- 
mark points and the bulk of the parameter space we consider (with DM masses in the 1 — 2 
TeV range). More generally, for DM masses in this range, regions where low- velocity capture 
is allowed are mostly already ruled out by constraints from the CMB, as shown in Figure 2. 

However, an scales very roughly as m x in models with the correct relic density, so 
the threshold mediator mass m t hrcs = a D m x/^ ^ m \- ^or m x ~ 2 TeV, if there are no 
additional annihilation channels to reduce the required value of ao, w t h rC s ^ 1 GeV, and 
it becomes critical to take WIMPonium formation into account for mediator masses in the 
entire sub-GeV range. 

For the purpose of studying WIMPonium formation, we restrict ourselves to the elastic 
case (5 = 0), as previously neglecting interactions involving the dark Higgs. As we have 
argued above, we expect the low-velocity capture cross section to trace the Sommerfeld 
enhancement, just with a different prefactor (higher by a factor of ~ 7.25). If this assumption 
is valid, the boost factor in the low-velocity limit still scales as anm x /m ( f ) oc m^/m^, so for 
close to the threshold value, BF sa t oc ttt," 1 , whereas the CMB constraint on this quantity 
scales as m x . Thus for mediators with masses calibrated to the WIMPonium threshold, the 
CMB limits become much less constraining at higher DM masses, and for m x > 2 TeV, 
we expect to find CMB-allowed regions of parameter space where WIMPonium formation 
cannot be neglected. 

Multi-TeV candidates to fit the cosmic-ray data must give rise to relatively soft e + e~ 
spectra with small or zero branching ratios to electrons, to avoid overproducing electrons 
at energies extending up to the DM mass; these features are not easily obtained from the 
simple vector portal models we have considered here (however, they can be accommodated in 
related models, e.g. [131]). In turn, such models require significantly higher present-day boost 
factors, as discussed earlier in the context of [34] - but conveniently, taking WIMPonium 
formation into account will inevitably give rise to just such larger boost factors. 



of s-channel annihilation through an off-shell 4>, hi analogy to true (ortho)muonium decay. However, in this 
analysis we will - as previously discussed for the 5 = case - neglect interactions between the <f> and other 
dark-sector states, such as the dark Higgs. 
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Figure 4. The present-day boost factor for the case of a single DM state with light enough to permit 
radiative capture into WIMPonium, as a function of the DM mass m x , and the mediator mass normalized 
to the WIMPonium binding energy. q_d is tuned to obtain the correct relic density. Red-hatched regions are 
ruled out at 95% confidence by constraints from WMAP5, taking the energy deposition fraction / = 0.2 to 
be conservative. 

For illustration, we solve the Boltzmann equation including WIMPonium formation 
and ^-dissociation for 2-4 TeV DM, assuming the WIMPonium capture rate to trace the 
Sommerfeld enhancement for v < ar> as discussed above, and taking the energy dependence 
of the (^-dissociation cross section to be [130]: 

4(1-7? cot- 1 »;) 1/n 

fJ ( £ ) Ke -i =2^T< V = (e-iy 1/2 , e = E/E thies . (A.6) 

1 — e Z7T " 

The results of this calculation are shown in Figure 4; we find that boost factors in the 600- 
700 range are possible if m x > 2.8 TeV, even for this elastic model (in comparison, the 
boost factor from Sommerfeld enhancement alone is ~ 100). Note that we neglect capture 
into excited states and transitions between the various bound states in the spectrum (for a 
discussion of such processes in the context of WIMPonium, see [132]); these processes are 
not generally negligible, but are beyond the scope of this study. In much of the parameter 
space for few- TeV DM, in order to evade the CMB limits the (ft mass must be quite close to 
the threshold value, and so capture into excited states will be kinematically forbidden at low 
velocities. 

B The Origin of the Mass Splitting in a Simple Model 

The high-energy Lagrangian in our "minimal model" takes the form, 

£ = (fy + igofa) ¥ + + ig D ^) h D (d„ - ig D <P») h* D - m x ^ 

~ 2T {** W D h h + ^°h D h D ) - \f»F% - e -F™F% + £ S M. (B.l) 
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We can use the high-energy limit to compute annihilation cross sections for the vari- 
ous processes that are approximately valid at all energies, since m^/m x is small and in all 
cases there are leading order terms with no dependence. The annihilation of fermions and 
antifermions into gauge bosons and oppositely-charged scalars have both been computed else- 
where and we will not give details here. The annihilation induced by the operator ^ c ^h* D h* D 
and its conjugate is easy to compute: the matrix element for e.g. fermion-fermion annihila- 
tion is just \M\ = \{y/K)u\V2\ (where "1" and "2" label the ingoing particles). Averaging 
over initial spins yields, 

(1/4) \M\ 2 = (l/4)(y/A) 2 Tr (pi^W - mim 2 ) = (y/A) 2 ( Pl • p 2 - mim 2 ), (B.2) 

S1,S2 

and dividing by two to account for the two identical particles in the final state, we find the 
COM frame cross section, 



2vr \M\ 2 1 s/2-2m 2 x (y ,2 \ , y 



° = ^-o / T = — - — ( -r ) Jl-Aml/s. (B.3) 



Now let us consider the symmetry breaking. We employ the notation of [133], and write 
^ as the Weyl fermion pair (%, rf). Then the ^ Cs $h* D h* D operator (and its conjugate) give 
rise to, 

Aspiit = (l/2){y/K)(xxh* D h* D + wh D h D + h.c). (B.4) 
Consequently, once /id develops a VEV, we can write the mass terms for x> V in the form, 

C n ^ = hx r,)(Z M ™*)( X )+h,., (B.5) 
2 \ m x m M J KV J 

where uim = (2//A)(/i_d) 2 . We can now write down the Takashi diagonalization matrix Q, and 
the mass eigenstates Xi>X2 : 

^=4f; \.y (?)=^( x ). (b.6) 



V2 V 1 ~V ' V X2 
Then the mass matrix becomes, 

2 V m X m M 



X2 



=k<* «>(" !> r M mv -,„„)(s)— ^ 

as desired, and we see that the mass splitting 5 = Imyi — 2(y/A)(/i£>) 2 . The splitting term 
in C transforms to, 

£ sp iit = (l/2)(y/A)((l/2)( X i + iX2?h* D h* D + (l/2)( X i - iX2?h D h D + h.c), 
= (l/4)(y/A)(xiXi(h D h D + h* D h* D ) - X2X2(h D h D + h* D h* D ) 
+ 2ixiX2{h D h D -h* D h* D )) + h.c. (B.8) 

Working in unitarity gauge, we will write ho — > {vd + p)/v2, so 5 = (y/A)v^; we then 
obtain, 

-Csplit -> (l/4)(i//A)(xiXi - X2X2)(v 2 D + p 2 + 2v DP ) + h.c, 
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so the strength of the XiXiPP interaction vertex is set by y/A (the Yukawa interaction is 
suppressed by vd/A). Furthermore, the mass of the gauge boson <fi is given by TnZ/2 = 
g 2 D {hr,) 2 , so = goVD- Thus we have the relation, 

y/A = 5/v 2 D = 6g 2 D /m% (B.9) 

Computing the rate for XiXi ~^ PP annihilation in the two-component formalism, we 
obtain, 

iM = 1 ( ^ ) ( V ^ & + X " X2a ) ' (B ' 10) 



2 



2 

I (pi-P2-my, (B.ll) 



Sg 2 



D_ 



exactly as previously. Since in the high-energy limit we can either say that half the particles 
are fermions and half antifermions, or that half are in the "ground state" and half in the 
"excited state" (really, the linear combinations of the Weyl fermions that will become the 
mass eigenstates at low energy), the — > h^liD and ipiipi — > m] cross sections should 
indeed be equal to obtain a consistent overall annihilation rate. 
Taking otp = gj^/An as usual, we can write, 

2 / x N 2 



a- * D 



4 1^2 V 1 " 4 ^/*- ( B - 12 ) 



In the low-energy limit s ~ Am 2 (1 + v 2 ), and v rc \ = 2v, so we obtain, 



™rel = ^ 2 — T ~ B - 13 

2 \ mi / mr. 
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